Utilization of the Wavefront Sensor and Short-Exposure Images 
for Simultaneous Estimation of Quasi-static Aberration and 

Exoplanet Intensity 

Richard A. Frazin 1 
rf raziri@umich . edu 

ABSTRACT 

Heretofore, the literature on exoplanet detection with corongraphic telescope 
systems has paid little attention to the information content of short exposures 
and methods of utilizing the measurements of adaptive optics wavefront sen- 
sors. This paper provides a framework for the incorporation of the wavefront 
sensor measurements in the context of observing modes in which the science 
camera takes millisecond exposures. In this formulation, the wavefront sensor 
measurements provide a means to jointly estimate the static speckle and the 
planetary signal. The ability to estimate planetary intensities in as little as few 
seconds has the potential to greatly improve the efficiency of exoplanet search 
surveys. For simplicity, the mathematical development assumes a simple optical 
system with an idealized Lyot coronagraph. Unlike currently used methods, in 
which increasing the observation time beyond a certain threshold is useless, this 
method produces estimates whose error covariances decrease more quickly than 
inversely proportional to the observation time. This is due to the fact that the 
estimates of the quasi-static aberrations are informed by a new random (but ap- 
proximately known) wavefront every millisecond. The method can be extended 
to include angular (due to diurnal field rotation) and spectral diversity. Numer- 
ical experiments are performed with wavefront data from the AEOS Adaptive 
Optics System sensing at 850 nm. These experiments assume a science camera 
wavelength A of 1.1 /i, that the measured wavefronts are exact, and a Gaussian 
approximation of shot-noise. The effects of detector read-out noise and other 
issues are left to future investigations. A number of static aberrations are intro- 
duced, including one with a spatial frequency exactly corresponding the planet 
location, which was at a distance of ~ 3X/D from the star. Using only 4 seconds 
of simulated observation time, a planetary intensity, of ~ 1 photon/ms, a stellar 
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intensity of ~ 10 5 photons/ms (contrast ratio 10 5 ), the short-exposure estimation 
method recovers the amplitudes static aberrations with a 1% accuracy, and the 
planet brightness with a with 20% accuracy. 

Subject headings: exoplanet detection 

1. Introduction and Motivation 

The idea of using millisecond exposures for ultra-high contrast imaging, such as is re- 
quired for exoplanet science, has been given little attention due to a variety of issues ranging 
from detector readout noise, to high-volume data management and undeveloped computa- 
tional methodologies. The method described here is to record the intensity in the science 
camera simultaneously with the wavefront sensor (WFS) measurements of the adaptive op- 
tics (AO) residual phase and solve an inverse problem to simultaneously determine the 
quasi-static (QS) aberrations and the planetary emission. Intuitively, it seems that the com- 
bination of millisecond exposures and the measurements of the WFS provide much more 
ability to estimate the QS aberrations and the planetary emission than do long-exposure 
measurements, even when combined with a scheme to estimate QS aberrations, such as a 
focal plane mask. There are two main reasons for this. The first reason is that the random 
AO residual phase provides a new pupil plane phase screen every millisecond, which com- 
pletely de-correlates on the time-scale of the wind blowing across the telescope aperture, the 
so-called "atmospheric clearing time," r c (Mcintosh et al. 2005). Thus, every r c seconds AO 
residual provides a statistically independent phase screen which modulates the quasi-static 
speckles in a new way. Each new phase screen provides more diversity in the observations 
and gives more power to estimate the QS aberrations and the spatial distribution of plan- 
etary emissions. In contrast, using a focal plane mask and with long-exposure images to 
estimate errors may introduce non-common path errors, and the focal plane mask provides 
only one measurement with no diversity. However, a focal plane mask could be combined 
with short-exposure imaging to increase the diversity, using a methodology similar to the 
one advocated here. 

The second reason for advocating millisecond images is that short-exposure measure- 
ments take advantage of the times when the random temporal modulation greatly reduces 
the intensity of the speckle background (at a given spatial point), which is the so-called 
so-called "dark speckle" principle described by Labeyrie (1995) and developed in the context 
of coronography by Boccaletti et al. (1998). The method described here can be considered a 
generalization of the dark-speckle method that takes advantage of both the WFS information 
and useful properties of the photon noise, particularly the fact that variance of the Poisson 
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distribution is equal to its mean (i.e., there is less shot-noise when the speckle background 
is faint). 

Roggemann and Meinhardt (1993) used wave- front sensor measurements to deconvolve 
images, but did not consider the effects of quasi-static speckle, which are critical for high- 
contrast imaging. Codona et al. (2008) describe a short-exposure technique called phase 
sorting interferometry which utilizes the WFS data to help determine the quasi-static aber- 
rations. However, unlike this paper, they do not consider planetary emissions, which have 
some potential to be confused with quasi-static aberrations. Gladysz et al. (2010) explored 
short-exposure statistics of stellar speckle and planetary emission, and showed that the 
resulting histograms provide some means to discriminate between planetary emission and 
speckle. However, reducing the intensity time-series to histograms discards much informa- 
tion, especially WFS' estimate of the AO residual phase, which is the function that is driving 
the temporal modulation in the first place. 

The current state-of-the-art post-processing methods, always applied to long exposure 
images, attempt to subtract the point spread function (PSF), which includes the effects of 
QS aberrations, by taking advantage of diversity introduced by the diurnal field rotation, 
multiple spectral channels and/or observations of naked stars with the telescope in a (hope- 
fully) similar state. The most well-known of these are called Angular Differential Imaging 
(ADI; Marois et al. 2006), and Locally Optimized Combination of Images (LOCI; Lafreniere 
et al. 2007). At the time of this writing, the latest variant along these lines utilizes principle 
components analysis to treat some of the difficulties of ADI and LOCI (including overly- 
aggressive subtraction and bias) and is called Karhunen Loeve Image Projection (KLIP) 
algorithm (Soummer et al. 2012). An independent line of thought in exoplanet detection 
involving the scanning Hotelling observer has been proposed by Caucci et al. (2007, 2009), 
who did not consider the effects of QS aberrations or a coronagraph (cf. Sec. [3] below). The 
work here is related to that of Sauvage et al. (2010), in which the authors proposed to use 
an analytical coronagraph model to simultaneously estimate the planetary emission and the 
stellar background from a long-exposure image. The most significant difference between this 
work and that of Sauvage et al. (2010) is the incorporation of WFS data and short-exposure 
images. 



2. A Simple Coronagraph Model 

Following Sauvage et al. (2010), consider a telescope with an AO system feeding a stellar 
coronagraph, shown schematically in Fig. [TJ Shown here is stellar variant of the classic Lyot 
coronagraph (Lyot 1939), versions of which have been used to study the Sun's corona for 
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many decades. In a Lyot coronagraph, light is first focused onto a focal-plane occulter, onto 
which the image of the star falls. The light that misses the occulter is then re-collimated 
and passed through the so-called Lyot stop, which removes much of the star light diffracted 
by the focal plane occulter. 

Assume an aperture function A(r), where r is the 2D spatial coordinate in the pupil 
plane. In this formulation, A(r) is assumed to include the effect of the Lyot stop in the 
coronagraph system (Sauvage et al. 2010) and may have a non-zero imaginary part. Above 
the atmosphere, the field due to the on-axis star is given by the (real) constant u*, and the 
field due to the planet is given by the linear phase-angle function u, exp[jka. p ■ r], where 
u, is real and u,/u* is the square-root of intensity contrast ratio between the star and the 
planet, k = 2ti/\ and A is the wavelength at the central wavelength of observation, and 
dp is the 2D direction vector (||aip|| < 1) indicating the location of the planet. This is 
easily generalized to multiple planets and extended objects (e.g., protoplanetary disks) with 
superposition of sources at angles {ck^}. The AO residual phase is given by <f) r (r,t) and the 
upstream (of the coronagraph) static aberrations (leading to so-called quasi-static speckle) 
are given by u (r). This formulation makes the simplification of excluding time-dependence 
of 4> u . While this is justified on the time-scales of seconds, or possibly minutes, it is not valid 
for longer observation times. One could describe slow changes of <p u , perhaps using B-splines 
to parameterize the temporal variation, but this is left to future investigations. 

The total wavefront phase aberration is given by: 



Generally speaking one expects that ||0 u (r)|| << ||0 r (r,£)||. The importance of M is due to 
the fact that it changes on very long time scales compared to <f) r and hence creates speckles 
that cannot be mitigated by long integrations. For simplicity, static aberrations downstream 
of the coronagraph, which are generally considered to be less important than upstream 
aberrations, are not considered, see Sauvage et al. (2010) for a discussion. The analysis 
here permits amplitude effects by allowing <p r . and <p u to be complex valued functions. This 
analysis assumes that the planet and the star are in close enough proximity so that they are 
in the same isoplananatic patch, but are separated enough for the planet to be beyond the 
influence of the coronagraph, i.e., 9 in < \\at p \\ < 9 iso , where 9 in is the inner- working-angle of 
the coronagraph and 9 iso is the isoplanatic angle. 

Without loss of generality, assume that the spatial average AO residual phase over the 
pupil is zero: 



as the value of r o does not have any effect on the intensity. Thus, the pre-coronagraph pupil 
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plane fields of the star and planet are given by: 

[r,t) = w*A(r)exp[j0(r,t)] 



u p (r,t) = u,A(r) exp\jkac p ■ r + j(f)(r, t)] . 



(3) 
(4) 



and the ~ superscript denotes that this field is upstream of the coronagraph (see Fig. [T]). 
Since the planetary light is assumed not to be effected by the coronagraph, this notation is 
omitted for u p . 

To model the effect of a the coronagraph in the optical system, one may use the for- 
malism of Sauvage et al. (2010) , in which "the perfect coronagraph is always defined as an 
optical device that subtracts a centered Airy pattern from the electromagnetic field." The 
perfect coronagraph concept seems to have originated with Cavarroc et al. (2006). In the 
pupil plane downstream of the coronagraph, the stellar field is given by 



u s (r,t) 



{T,t)-U*T){t)A{T) 



u* {A(r) exp[j'0(r, t)] - ri(t)A(r)} 



(5) 
(6) 



where r)(t) is the complex number (|| rj(t) ||< 1) that minimizes the total energy transmitted 
by the coronagraph, J u s (r, t)u*(r, t)d 2 r . It is easy to show that 

/ d 2 rA(r)A*(r)exp[-j0 r (r,t)] 



7}(t) 



(7) 



/ d 2 rA(r)A*(r) 

where the * superscript indicates complex conjugation, and the u (r) has been ignored as it 
has little consequence here. Note that as 0(r,t) — > 0, rj(t) — > 1, implying the extinction of 
all the stellar emission, which is the defining characteristic of a perfect coronagraph. 

Under Fraunhofer diffraction, which relates image plane fields to pupil plane fields by 
Fourier transformation (e.g., Goodman 1968), the pupil-plane fields u p and u s lead to the 
following post-coronagraph, image-plane fields for the planet and star: 



v p (p, t) 
v a (p, t) 



. J d 2 r A(r) exp j 
d 2 r A(r) exp 



-fp- 

k 

- }J pr 



fa. p ) ■ r + r (r,t) 
exp[?'0(r, t)] -rj(t) 



(9) 



where p is the 2D image plane coordinate, / is the effective focal length, and 4> u (r) has 
been dropped from Eq. (JS}, as the planetary light scattered by the static aberrations is not 
important. The instantaneous intensity in the image plane is given by the incoherent sum 
of the instantaneous planetary and stellar intensities: 



I(p,t) = I s {p,t) + Up,t) 



(10) 
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where I p (p,t) = v p (p,t)v*(p,t) and I s (p,t) = v s (p,t)v*(p,t). Using Eq. (JSJ), the planetary 
intensity is given by: 



uz 



expj 



d 2 r / dVA(r)A*(r') x 



k 



-(p-fcx p ).(r'-r) + ( f> r (r,t)-4>* r (r',t) 



11) 



Similarly, from Eq. ([9]), the star's intensity is given by: 



I S (P,t) 



u: 



d 2 r / dVA(r)A*(r')exp 



x 



77(*)»7*(*) + e J ' [ * (r ' 4) -^ (r ''* )] - r/(t)e-^ (r ''' ) - ^(t)e^ (r '' } } , (12) 



where the reader is reminded that (j)(r,t) = <j> r (r,t) + u (r). From this equation, it can 
be seen that the star's instantaneous intensity is a complicated function of AO residual 
<j) r (r,t), the static aberrations 4> u (r), and the complex extinction factor rj(t), which accounts 
for the effect of the ideal coronagraph. In contrast, the planetary intensity is much more 
simple because the static aberrations were dropped and it was assumed earlier that the effect 
of the coronagraph is not important since ct p is greater than the inner working angle 9i n . 
Sauvage et al. (2010) obtained an expression equivalent to Eq. (1X21) . except their analysis 
also included aberrations downstream of the coronagraph. As Sauvage et al. (2010) were 
primarily concerned with long exposure images, they took the time average of Eq. ()12p . 
which eliminates (j) r (r,t) in favor of structure functions of the AO residual. 



2.1. multiple planets and/or extended emission 



Equation (11 II) assumes planetary emission only from the sky-angle a. p , but it easily 
generalizes to the case of extended planetary emission and multiple planets. To account for 
such emission, one must sum incoherently over all of the various sources. Let J (a) be the 
brightness distribution on the sky (excluding the star), then the generalization of Eq. ( TTTT) 
is: 



Ip{p,t) 



d'aJ(a) / d"r / dVi(r)i*(r') x 



expj 



p-f a ).(r'-r) + <t> r (r,t)-<t>*(r',t) 



(13) 



Henceforth, the double integral over r and r' in Eq. (JTTJ will be called the " planetary 
intensity kernel." 
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2.2. modulation of planetary intensity and stellar speckle 



The temporal variation of the AO residual modulates the stellar speckle and the plane- 
tary emission differently and in a way that can be exploited by short-exposure imaging. At 
the planet location in the image plane, p = fct p , consider a 1st order Taylor expansion of 
the exponential in Eq. flH]), which gives: 



where the equality makes use of Eq. ([21 . Thus, the intensity has no terms that are linear in 
4> r , and, to first order, the planetary field isn't modulated by the temporal variation of the 
AO residual phase. Temporal variability in the planetary intensity is due to the higher order 
terms in the expansion. In contrast, at the point p = fct p , Taylor expansion of the second 
exponential in Eq. (J5J) yields: 



One can see that Eq. ( !T5|) indicates that the field should be greatly modulated by the residual 
atmospheric fluctuations at the location of the planet. The key differences in the temporal 
behavior between the planetary case [Eq. ( fl4]) ]. and the stellar case [Eq. ( |T5l) ] is due to 
the spatially oscillating exp(— jkct p ■ r) factor inside the integrand, and the fact that 77 is 
a function of time. Fig. [2] gives an example, taken from Experiment 1 in Sec. HJ of the 
difference in the modulation of the planetary vs. stellar brightness at the same point in 
the image plane. In this case, the stellar brightness is enhanced by a static aberration with 
the spatial frequency corresponding to the planetary position. Figure |3] shows histograms 
of the planetary intensity (top) and the stellar intensity (bottom) from the same time- 
series. The positive skewness of stellar distribution and the negative one of the planetary 
distribution was noted by Gladysz et al. (2010), who proposed using these properties to 
separate planetary light from the stellar speckle. In the case shown in Figure El the stellar 
intensity was a nearly perfect negative exponential distribution, corresponding to so-called 
"fully developed speckle" (Goodman 2007), but other examples not shown here indicate that 
the stellar intensity may deviate from this distribution substantially. 

The differences between the planetary and stellar modulation implied by Eqs. ([H]) and 
09]) or Eqs. ( }T4"j) and ffl~5]) hint that short exposures may provide a means to discriminate 
between between the planetary emission and speckle from the star. Thus, one would expect 
that the WFS estimate of </> r (r,i) should provide some additional means to separate the 
planetary emission from the speckle. 




(14) 



v s (fcx p , t) « / d 2 r A(r) exp [jka p ■ r] 1 



(15) 
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2.3. Series expansion of the QS aberrations 

Our problem requires estimation of the static aberrations, which is less difficult when 
they can be assumed small enough for the following linearization to be valid: 

exp[70 u (r)] « 1 + j<f> u (r) . (16) 

This is not particularly stringent, as one expects that static errors should satisfy ||^» u (r)|| << 
1 in a good optical system. Note that the AO residual is not assumed to be small and eP^ T,t ^ 
is not linearized in this paper. Calculating the instantaneous total intensity from Eq. (fTUl) 
using Eqs. ©, flU, and flU} yields: 

I(p,t) = I p (p,t) + A(p,t) + (Bo</>*)(p,t) + 

(B*o0 tt )(p,t) + (Co[0 tt ,^])(p,t) (17) 

where the o symbol indicates composition. Of course, this equation may be generalized to 
include extended objects and/or multiple planets by replacing the first term with Eq. (fl3"|) . 
The linearization in Eq. (|T6|) applied to Eq. (1T2"]) . creates three types of terms: those that 
do not depend on </> u , those that are linear in <p u and those that are quadratic in <f) u . These 
terms are represented, respectively, by the A, B, and C terms in Eq. (fTTj) . The B terms lead 
to the so-called "pinned speckle" in long exposure imaging (Bloemhof 2004). Before giving 
formulae for these operators it is helpful to define the kernel function: 

.k , 



G(r,r'; p) = A{t)A*{t') exp 
which has the useful symmetry property: G(r, r'; p) = G*(r', r; p). 



A(p,t) = ulJdhjdh'G(r,r' ] p){r ] (t)i 1 %t)+expjiMr,t)-4 > *(r',t)] 

-r)(t)exp[-j(t>*(r',t)} - ^(t) exp[j'0 r (r,t)]} (19) 
(Bo^)(p,*) = ulj d 2 r j dVG*(r,r';p)f u (r) x 

{jr](t)exp[-j<f>*(r,t)} - j exp j[0 r (r', t) -0*(r,t)]} (20) 
(^o0 u )(p,t) = ulj d 2 r y dVG(r,r';p)0 u (r) x 

{ - j V *(t) exp[j0 r (r, t)] + j exp j [0 r (r, t) - 0^(r', t)] } (21) 
(Co[0 M ,0^])(p,t) = ^1 d 2 r y dVG(r,r';p)0 u (r)0:(r')expj[0 r (r,t)-0*(r',t)](22) 
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In Eq. (j!7p . the A term is due to the atmospheric residual speckles, the C term is due to the 
static aberrations (which are also modulated by the AO residual), and the B terms are due 
to the mixing of these two effects. Note that when <p u = 0, I s {p, t) = A(p, t). While all of the 
integrals in Eqs. (fT9|) through f ]2"2"]) can be separated to integrate over r and r' independently, 
the form given here is more conducive to presentation. 

To make estimation of the QS aberrations tractable, a finite expansion for (f) u (r) is 
assumed: 

K 

u (r) = ^a fe ^(r) (23) 
fe=i 

where K is the number of terms required and {^fc( r )} are basis functions on which the 
static aberrations can be represented. Examples of possible expansion bases include Fourier, 
Zernicke, and B-splines (e.g., Unser 1999), which include pixel bases. Importantly, orthogo- 
nality of the {ipk( r )} w iH n °t be assumed, allowing considerable freedom in the form of the 
expansion. For convenience, the quantities {dk} ar e placed into a 1 x K (column) vector a. 
Substituting Eq. (123]) into Eq. ( 1T7|) gives three types of terms, those that are independent of 
a, linear in a, and quadratic in a: 

I(p,t) = u 2 ,i p (p,t) + A(p,t) + 

a H b(p, t) + b H (p, t)a + a H C(p, t)a , (24) 

where the H superscript indicates the transpose and complex conjugate (Hermitian conju- 
gate), i p (p,t) = I p (p,t)/ul is the planetary PSF. Eq. ( 124"]) expresses the instantaneous inten- 
sity in terms of the unknown planetary intensity u 2 , and the static aberration coefficients a. 
The column vector b(p, t) has K components, and C(p, t) is a Hermitian symmetric, K x K 
matrix. Thus, for fixed (p,t), a H b(p, t) is the scalar product of the two complex vectors. 
Similarly, C(p, t)a is the standard matrix-vector multiplication, with a H C(p, t)a resulting 
in a scalar. Using Eqs. (120]) . (|22|) and (|23|) . the components of b and C are given by: 

[b(p,t)] k = u\J d 2 r J dVG*(r,r';p)^(r) x 

{jr)(t)ex V [-j(f>* r {r,t)} - jexp j[0 r (r',t) -<#(r,t)]} (25) 

[C(p,f)]« = ulj d 2 r J dVG(ry ;P )^(r)^*(rO x 

expj^^t)-^',*)}. (26) 

Note that Eqs. fT25l) and (126]) are separable and can be computed rapidly with fast Fourier 
transforms (FFTs), as can A(p,t). 
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3. Detection/Estimation Problem 

The problem of detecting and estimating the planetary emission in the framework pre- 
sented here can be posed as separating the planetary and stellar emission given noisy mea- 
surements of I(p,t) and the WFS data, summarized as the estimate of the AO residual 
4> r (r,t). It follows that using the WFS to estimate <f) r (r,t) is a key aspect of the proposed 
methodology^ While some current WFS' estimate other quantities, such as V0 r (r, £) in the 
case of Shack-Hartmann sensors, estimating <p r (r,t) from the WFS data is always possible 
with additional processing. Optimal estimation of <p r (p,t) with the maximum likelihood 
method is discussed in (Barrett et al. 2007) and (Bechet et al. 2009). 

The most general case allows the planetary emission to be a function of the sky-angle a, 
and then the problem would be estimate the spatial function J(ot) simultaneously with the 
QS aberration coefficients a, a problem which falls under the category of estimation theory. 
In contrast, the problem of trying to decide whether or not there is planetary emission falls 
under the regime of detection theory. Considering the problem of exoplanet detection within 
the framework of classical detection theory, Caucci et al. (2007) simulated the scanning 
Hotelling observer (SHO). The SHO is an example of a "scanning detector," in which one 
considers the probability of detection of a source at one location otk at a time. The process 
of evaluating the a series of locations is called "scanning." In Caucci et al. (2009), they 
generalized their previous analysis to include a time-series of images. The possibility of 
multiple planets considerably complicates scanning observer framework and it is not suited 
for estimation of extended emission. The, perhaps naive, first step approach in this paper 
is to assume that emissions other located outside of the scanning template, i.e., the group 
of pixels simultaneously searched for planetary emission, do not significantly influence the 
estimate of the planetary emission or its error bar. The numerical examples in Sec. H] support 
this notion. With this approach, planetary emission that is estimated to be non-zero with 
sufficient statistical confidence can be considered a detection. 



3.1. Linear System Formulation 

Here, Eq. (1241) is put into the canonical form for linear systems, y = Hx, where y 
represents the "observation" vector, x is the vector of unknowns and H is the system model 
or "data," as it is sometimes called in the statistics literature. 



1 Alternatively, one may also consider ignoring the WFS' estimate of 4> r (r,t), and instead, find suitable 
estimators for A, B, and C directly from the WFS data. 
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While Eq. fl2l]) has only K + 1, unknowns, it is nonlinear due to the quadratic term, 
a H C(p, t)a, which may introduce numerical difficulties. In exchange for creating an addi- 
tional set of unknowns, which will increase the uncertainty in the resulting estimates, one 
may remove this nonlinearity. The new K 2 unknowns are the quanties {a^al}, which are 
placed into the 1 x K 2 vector d. Similarly rearranging the matrix elements C(p, t) into the 
K 2 x 1 vector c H (p, t), one obtains: 

a H C(p,t)a = c H (p,t)d. (27) 



Consider a spatial template with iV locations {pi, Pn}, and take pi = fa p to be the 
location at which one desires to know whether or not a planet is present and how bright it 
might be. It is assumed that the science camera exposures are simultaneous with the output 
of the AO system wavefront measurements, with T exposures taken at times {ti, tr}- 
These exposure times form the temporal template. Thus, there are NT observations to be 
placed into the vector: 

r yi 



(28) 



where the {y,} vectors (one for each exposure) are given by: 



" I(pl,ti) 




' A{pi,U) ' 






_ A{p N ,ti) _ 



(29) 



Similarly, the (column) vector x = [u 2 , a T , a H , d T ] T , where T indicates transposition (without 
complex conjugation). Following the same notation scheme, 



H 



(30) 



where the individual matrices {Hj} have the structure: 

i p (pi,U) b T (pi,ti) b H (pi,tj) c H (pi,tj 



i P (pN,ti) b T (p N ,U) b u (p N ,U) c H (pi,ti 



(31) 
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Then, under the linear formulation, the estimation problem is to find the values oi ul, a and 
d which best satisfy: 

(32) 

In the classical form for linear systems of equations, Eq. (132 p expresses the relationship be- 
tween measured instantaneous intensities, which are one of the terms in y, and the unknown 
planetary intensity ul, the static aberration coefficients a, and the products thereof, con- 
tained in d. Note that the scale of the linear computation may not be especially daunting, 
as there are only K 2 + 2K + 1 unknowns, which are much less than the NT data points. For 
example, the matrix H H H only has K 2 + 2K + 1 singular values. 

When considering spatially extended planetary emissions, the x vector is expanded 
so that ul is replaced a column vector representing the extended emissions and i p (pi,t) 
is replaced by a matching row vector accounting for the convolution with the planetary 
intensity kernel [cf. Eq. ( ITS]) ]. 

If a is complex, then the estimate can be improved by demanding consistency between 
the estimates of a and a*, or, re-expressing the system in terms of real and imaginary parts. 
Of course, if a is assumed to be real, meaning static errors only affect the phase, not the 
amplitude, x simplifies to x = [ul, a T , d T ] T , and the 2 nd and 3rd columns of H are summed 
together to form a single column. In addition, a positivity constraint may be applied to the 
planetary intensity. 

3.2. Quadratic System Formulation 

As explained at the beginning of Sec. 13.11 the linear system in Eq. fl3"2"|) requires K 2 + K 
aberration coefficients instead of only K, and both formulations have the same input data 
from which these parameters must be estimated. It follows that the quadratic formulation 
will produce more accurate estimates with tighter error bars, simply because the number 
of quantities that must jointly estimated is much smaller. The only reason not to employ 
the quadratic formulation is due to the numerical difficulties that may be encountered when 
solving the system and estimating the uncertainties. The linear formulation may improve 
the computational possibilities for the quadratic system by providing a "warm start" (i.e., a 
good initial guess). It may be useful to note that given linear estimates of a and d, a and d, 
and their error covariances C aa , C a d and Cdd, in a second computation, one could improve 
the estimates of a by using the knowledge that each component of d is equal to a«a* for a 
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particular choice of i and j. These improved estimates may provide a better warm start for 
the quadratic system solver. 

The objective of remainder of this section is to write the system in the form y = 
H'x' + Q(x'), where Q is the quadratic term. The vector y is the same in both the linear 
and quadratic formulations. The quadratic version of H, H' is the same as H, except with 
the final column removed: 



H' 



i p (pi,U) b T (pi,ti) b H (pi,tj) 
i P (pN,ti) b T (p N ,ti) h H (p N ,ti) 



(33) 



so that H' = [Hj, H£] T Similarly, x' = [u 2 „ a T , a Hl 



To put the quadratic term in a form suitable for digital computation, first note that for 
each of the NT points in the spatio-temporal template one must calculate the scalar value 
of a H C(p, t)a, which is achieved with the following formulation: 



Q(a) = (ljv T ®a H ) 



C{p u h) 

C(p N ,ti) 
C(pi,t 2 ) 

_ C(p N ,t T ) _ 



(34) 



where Itvt is the row vector [1, 1] with NT components, and ® is the Kronecker matrix 
product. This Kronecker product creates a large row vector with NT copies of the row vector 
a H . Then, the quadratic variant of the problem can be stated as finding the u 2 m and a which 
best satisfy: 



H' 



a 

a* 



(l WT (8)a H ) 



C(pjv,*i) 
C(pi,t 2 ) 

_ C(p N ,t T ) _ 



(35) 



Eq. (l35il is the quadratic analog of Eq.(l32l). 
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3.3. Uncertainty Propagation 

This section considers uncertainty propagation within the linear formulation of the prob- 
lem, y = Hx provided in Sec. 13.11 The quadratic version of the problem, formulated in 
Sec. I3.2[ will have smaller uncertainties due to the much smaller number of variables to 
estimate, but determining them requires more specialized mathematical procedures than the 
linear case. 

This section will assume that NT is much greater than the number of parameters to 
be estimated. In practice this will be satisfied, as millisecond exposure times will produce 
a lot of data. For example, with a spatial template consisting of 1000 pixels, and exposure 
times of 10 -3 seconds, one hour of observation gives NT = 3.6 x 10 9 data points. If the 
QS aberration expansion has, say, 20 terms, the x vector will have 441 components, and the 
linear problem can be considered with the framework of least-squares. 

The first term in y, the observed intensity /, inherits its statistical properties from 
noise in the observations (most likely dominated by Poisson counting noise and detector 
read noise), which are treated with relative ease. However, it is of critical importance to 
bear in mind that the errors in (j) r (r,t) will cause uncertainties in both y and H which 
must be propagated to x [cf. Eqs. (|2"B|) and (I3"0"j) ]. Standard least squares is likely to be 
insufficient, and total least squares (TLS) is a more appropriate analysis framework (Van 
Huffel and Vanderwalle 1991). 

The short-exposure approach proposed here is unique among the various processing 
methods in that knowledge of the quasi-static aberrations increases with the observation 
time, whereas other methods are limited by the systematic error of imperfect a-priori knowl- 
edge of the PSF. As the condition of the matrix H improves with each random wavefront, 
increasing observation time reduces the covariance of x more rapidly than 1/T. This can 
be shown for the simple case of a perfect WFS, for which (j) r (r,t) = r (r,t). This implies 
that H is known exactly, as is A(p, t), so that the only uncertainty remaining comes from 
the measurement noise in I(p,t), which is manifested in y according to Eq. (129 p . Consider 
a spatial template of N detector pixels and assume that the measurement noise of each 
exposure has an N x N diagonal covariance matrix Ejv, which is the same for each of the T 
exposures. Then, the total covariance matrix for all of the NT data points is given by large 
diagonal matrix S = It eg) Sjy, where It is the T-by-T identity matrix, and log-likelihood 
function is therefore given by: 

lnL(x|y) ex -(y - Hx) H S" 1 (y - Hx) . (36) 

This likelihood is maximized by x = (H H S _1 H) _1 H H S _1 y, and the covariance of x is given 
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by the (K 2 + 2K + l)-by-(K 2 + 2K+1) matrix: 

E x =(H H E- 1 H)- 1 =ff;H?E^H i ) . (37) 



i=i 



Now, each of the matrices (Hf E^- Hi) is positive semi-definite, and the trace (which is 
the sum eigenvalues) of the sum ^^ =1 (H?E^ 1 H i ) must increase, on average, linearly with 
the number of exposures T. Furthermore, since each of the Hj is a manifestation of an 
independent, random phase screen <f> r (r, ti), the matrix [H^, H^ r ) _ 1 ] T , will have an equal or, 
more likely, smaller ratio of largest-to-smallest singular values than the matrix [ttj, H^] T 
(in other words, the condition of H improves with increasing observation time). Thus, one 
has: 

S * = (X>^ lR *) < ( T < H f S ^ H 0) _1 > (38) 

where the brackets < > indicates ensemble average over the statistics of the AO residual. 
Thus, it is clear that E x decreases more quickly than 1/T. This better than 1/T perfor- 
mance was observed in the numerical experiments in Section HI which compared simulated 
observations of 1 and 4 s. Demonstrating the inequality Equation fl38l) in a rigorous manner 
is beyond the scope of this paper, but the mathematically inclined reader may consult Tropp 
(2012). While generalization of this proof to the case in which H also has uncertainty and 
its own covariance matrix is likely to be substantially more difficult, it seems that this basic 
result is unlikely to change. 



4. Numerical Simulations 



This section explores some of the concepts described above with numerical simulation. 
In the experiments described below, the residual phase data taken by the WFS on the 
AEOS Adaptive Optics System (Roberts & Neyman 2002), interpolated onto an 80 x 80 
grid with a circular aperture inscribed (the radius D/2 of the circle corresponding to the 
telescope aperture was 40 pixels), was used as r (r,t). The interpolation was done with a 
standard spline interpolation routine. In this way, realistic wavefronts were used. By taking 
the WFS data as the true wavefront, the problem imperfect knowledge of the wavefronts 
was eliminated in these proof-of-principle experiments. Including the effects of imperfect 
wavefront knowledge is the subject of a forthcoming study. AEOS has 941 actuators and the 
wavefront data utilized here corresponds to a sampling rate of 10 -3 s. The AEOS sensing 
wavelength was 850 nm and these numerical exercises assumed a science camera operating at 
1.1 fi that samples at exactly the rate as the WFS (with no delay). The residual phase value 
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was converted from the 850 nm value to the 1.1 /i value by assuming that the optical path 
length difference is independent of wavelength, so that the residual phase in radians simply 
scales inversely with the wavelength. In reality, the 1.1 fi wavefront would have amplitude 
errors and non-common path errors, but such effects were not considered. 

As per Equations ([3]) and (J4]), the planetary and stellar pupil plane fields were subjected 
to the same static aberrations as well as the random AO residual. The static aberration was 
also placed onto the 80 x 80 grid with a circular aperture inscribed. The static aberrations 
included 11 terms in the expansion shown in Eq. (123 p : Zernike functions numbered 41-50 as 
well as a sinusoid which places a speckle at the planetary location as well as its conjugate 
point. The coefficients {a^} as well as the aberration functions {ipk} were real, corresponding 
to phase-only aberrations. The chosen values of {a^} resulted in the aberration function 
shown in Figure [5J All Fourier transforms were calculated by zero-padding the 80 x 80 arrays 
into a 256 x 256 arrays and then using the FFT. Figure |6] shows the central portion of the 
so-calculated magnitude of the Fourier transform of static aberrations shown in Figure |5j 
The two bright points are the result of the sinusoidal aberration, and the bottom one is 
exactly spatially coincident with the planet location in Experiments 1 and 2 below. With the 
aperture and zero-padding as described above, the minimum of the Airy function 1.22X/D 
when calculated analytically) is about 4 pixels, so, in these simulations, 1 pixel (in the image 
plane) corresponds to an angle of about 0.3 X/D. The bottom spot in Figure [5] and the 
planet are 10 pixels away from the center of the image (the star's location), corresponding 
to a distance of ~ 3X/D. 

The results of 3 experiments, numbered 1, 2 and 3, are given below. The purpose of 
Experiment 1 is to demonstrate the ability to detect a planet under the above conditions with 
an assumed position (pixel location). Treating the more common case of an unknown planet 
position would require a combination of generalizing the formulation to consider emission in 
a group pixels, called the "scanning template" , as described near the end of Sec. 13. 1\ and, 
in a series of subsequent calculations, moving the scanning template to cover the required 
region, as described near the beginning of Sec. |3j The experiments here correspond to a 
scanning template consisting of a single pixel. Experiment 2 identical is to experiment 1, 
except that it uses only one-quarter of the exposures in Experiment 1. The purpose is, upon 
comparison to Experiment 1, to provide a demonstration of the claim in Equation [311 which 
is that the estimate variance decreases faster than than 1 /T, where T is the total observation 
time. Experiment 3 concerns the viability of the scanning concept. Ideally, emission from 
outside the scanning template would not influence the estimation of emission from within 
the template, but, in practice, there will be leakage. This experiment gives an indication 
that leakage is unlikely to be a severe problem. All experiments used the linear formulation 
in Section 13.11 Presumably, utilizing the quadratic formulation of the problem described 
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in Sec. 13.21 would lead to superior results (due to the smaller number of unknowns) than 
those reported here, but the difficulties in minimizing the quadratic form in Eq. (135]) and in 
determining the uncertainties in the estimates are left to a later study. 

These experiments utilized either 1 s (1000 exposures) or 4 s (4000 exposures) of AEOS 
data. Figure H] shows the 1.1 /x Strehl ratio for all 4000 exposures, given by \rj(t)\ in Eq. 
Also depicted on the same graph is the fraction of the stellar emission transmitted by the 
coronagraph. In the very best frames, it removes about 95% of the stellar emission and much 
less when the Strehl ratio falls. The 850 nm Strehl ratios are smaller and are not shown 
here. The contrast ratio between the planetary and stellar brightness was 10~ 5 , and the 
stellar intensity about 10 5 photons/ms. Integrated over the entire image plane, the time- 
average planetary intensity was about 1 photon/ms, however, at the pixel corresponding to 
the center of the planet's Airy pattern, the intensity was only about 0.045 photons/ms/pixel. 
The average stellar intensity at that pixel (enhanced by the sinusoidal aberration) was over 
500 times greater, with an intensity of over 23 photons/ms/pixel. Integrated over the image 
plane, the time-average stellar intensity was about 4.5 x 10 4 photons/ms, meaning that, 
on the average, the coronagraph killed about 55% of the star light. The noise-free stellar 
brightness, averaged over all 4000 frames, is shown Figure where planet is too faint to be 
seen. 

The numerical experiments here simulate the shot-noise that is associated with photon 
detection, but do not include detector noise, which is left to a forthcoming study. For the 
purposes of the signal modeling in this paper, the semi-classical model for photoelectric 
detection (Goodman 1985) is adequate. In the semi-classical photo-detection paradigm, 
there is a classical intensity (i.e., not quantized) incident upon the detector, but the detector 
measurement is modeled as a Poisson distributed random variable with the mean equal to the 
classical intensityl In order to avoid the complications associated with estimation in non- 
Gaussian noise, Gaussian noise was chosen to behave similarly to the Poisson distribution. 
Thus, if the (noise-free) intensity was I photons in a millisecond exposure, the "measured" 
value (i.e., the simulated measurement) would be I + n, where n was taken from a zero-mean 
normal distribution with a variance equal to /. This Gaussian shot-noise was added to each 
pixel in each frame, constituting a single noise realization. Additional noise realizations 
within the same experiment used same intensity values /, but different random values of n. 
The experiments were run multiple times, with different noise realizations in order to compare 
the variances in the estimates with the analytically derived variances from Equation (I37p . 

As per Equation (|2"9~I) . a spatial template {p 1 , p^} was chosen. For these experiments, 



2 The units of the classical intensity must be the [photons/time] x [exposure time] . 
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this template corresponded to the central 27 x 27 pixels of the 256 x 256 zero-padded array 
in which the Fourier transforms were calculated. This template choice excluded 30% of the 
stellar photons and 20% of the planetary photons from analysis (this difference is probably 
due to the coronagraph) , resulting in an H matrix (for the 4000 exposure case) with 4000 x 
27 2 = 2916000 rows and 1 + 11 + ll 2 = 133 columns. 

Experiment 1 included 4000 exposures, with a mean Strehl ratio of 0.75, and was run 
10 times. The only difference between the runs was the realization of the noise, with each 
noise realization consisting of 2916000 random numbers. These 10 runs had a mean estimate 
of the planetary intensity of 1.1 photons/ms and a standard deviation of 0.16, which is 
comparable to the analytical result [i.e., the square-root of the 1-1 element of the estimate 
error covariance matrix in Equation 037p ] of 0.2. Fig. [S] shows the estimates from one of 
the runs and the corresponding true values of the linear coefficients {a^}- The error bars, 
which are smaller than the size of data points, come from the square-root of the diagonal 
of the estimator covariance matrix. Similarly, Fig. [9] shows the estimates and true values 
of the quadratic coefficients {a^a*} (although complex conjugation was not needed since all 
coefficients were real). 

Experiment 2 was the same as experiment 1, except with only 1000 exposures, with a 
mean Strehl ratio of 0.76, were used. Since Experiment 2 used one-quarter of the exposures 
used in Experiment 1, a simple 1/T dependence of the estimator variance would lead to 
standard deviations that are double those of Experiment 1. Experiment 2 was run 20 times, 
each with a different realization of the random noise described above. The mean estimate of 
the planetary intensity in these runs was 0.7 photons/ms and a standard deviation was 0.96, 
comparable to the analytical result of 1.1. Comparing the analytical standard deviations 
from Experiments 1 and 2, one can see that 1.1 >> 2 x 0.2, supporting the contention 
in Equation (1381) that the variance of the estimate decreases more quickly than 1/T. The 
estimates of the static aberration coefficients had analytical standard deviations about 6 
times those in Experiment 1. The equivalent to Fig. [8] is visually undistinguishable, as 
the error bars are still quite small. The average analytical variance of the 121 quadratic 
aberration coefficients was about 2.7 times that of those in Experiment 1, again supporting 
Equation fl3"g|) . 

Experiment 3 was designed to test the algorithm's stability to planetary emission from 
angles other than cx p , to help evaluate its utility as a scanning detector, in which one moves 
the scanning template and reprocesses the data for each template location. Ideally, planetary 
emission from outside the scanning template would have a negligible effect on the estimates 
of both the static aberration coefficients and the planetary emission at cx p . Experiment 
3 was the same as Experiment 2, except the planet was placed at an angular distance of 
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about 1.75X/D (7 pixels on the image grid) to the right of the assumed planetary position 
(i.e., that of Experiments 1 and 2). The same matrix H from Experiment 2 was used, 
which assumed any planetary emission was coming from the planet location in Experiments 
1 and 2 (see Figure Thus, the matrix H assumed that the planet was in the wrong 
location. In addition, the planet's brightness was increased to 1000 photons/ms, reducing 
the contrast ratio to 10 2 , and making the planet easily visible in the time-averaged image. 
The ideal result of this experiment would be a planetary intensity estimate of and the 
correct estimates of the static aberration coefficients. The estimated planetary intensity 
was about 25 photons/ms, or 2.5% of the true planetary brightness. Thus, this experiment 
showed 2.5% leakage from a distance of 1.75X/D, and the scanning detector response is 
strongly peaked at the correct location. Importantly, in this experiment, the estimates of the 
aberration coefficients were nearly as good as in Experiment 2. When the planet was placed 
at the smaller distance of about 0.9X/D, the leakage was much greater - about 12% and the 
static aberration coefficients where not as accurate. This indicates that scanning template 
of a size of several X/D should be sufficient to implement a scanning detector/estimator. 
Such a scanning detector could be run iteratively, with the second iteration using a spatial 
template that includes all likely source locations from the first. 



5. Conclusions 



The purpose of this paper is to illustrate the possible potential for short-exposure imag- 
ing of exoplanets and thereby spur further research on the topic of short exposure imaging 
for high contrast applications. Specifically, it describes a method for simultaneously esti- 
mating the static aberrations and the planetary emission from millisecond exposures with 
ground-based stellar coronagraphs. The ability to estimate planetary intensities in as lit- 
tle as few seconds has the potential to greatly improve the efficiency of exoplanet search 
surveys. The new and crucial aspect of this work is that the method incorporates the wave- 
front sensor's estimate of the AO residual phase to provide much needed diversity in the 
data. Thus, new and statistically independent information about the planetary emissions 
and quasi-static aberrations is revealed on the atmospheric clearing time-scale (Mcintosh et 
al. 2005), which is, more-or-less, the telescope diameter divided by the wind speed. It is this 
diversity property that allows the covariance of the estimate of the planetary emission and 
the QS aberrations to decrease more quickly than inversely proportionally to the observation 
time. While the random AO residual is an extremely useful source of diversity that can be 
used to characterize the optical system while simultaneously estimating the planetary emis- 
sion, a fundamental trade-off is shown in Fig. HJ Large AO residual is accompanied by poor 
rejection of stellar emission by the coronagraph. 
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Due to the exploratory and preliminary nature of this work, many assumptions were 
made for the sake of simplicity. Among these assumptions are the the lack of aberrations 
downstream of the coronagraph, an ideal coronagraph, Fraunhofer diffraction, and no tem- 
poral variation in the static aberrations, to name a few. The model can be generalized to 
incorporate more effects. For example, it can allow for temporal variation of the quasi-static 
aberrations on longer time-scales using a scheme such spline interpolation in the time dimen- 
sion, at the expense of having to estimate the spline coefficients. Furthermore, the method 
given here generalizes to the include spectrally resolved data and the diurnal field rotation. 
E.g., to take advantage of S spectral channels, one must solve for S planetary intensities 
simultaneously, but there are still only K independent values characterizing the QS aber- 
rations, since the phase variation scales with wavelength. Similarly, to take advantage of 
diurnal field rotation (angular diversity), one only has to apply a coordinate transformation 
before including the additional data. 

The fundamental limits of the short exposure approach are not easy to determine, as 
they are set by the accuracy of the various modeling aspects of the coronagraph system, 
including the reduction of the QS aberrations to a finite set of parameters, as well as the 
various/bias characteristics of the estimator of the AO residual. In a practical application, 
these will need to be evaluated with care. 

In order to apply the method proposed here to real data sets a number of challenges 
must be met. While any practical realization will no doubt have to consider various details, 
e.g., time differences in the exposures of the WFS and the science camera, the more essential 
tasks include the following: 

• Develop detectors with low readout noise so that short exposures are not penalized 
heavily. Shorter wavelengths at which readout noise is less of a challenge should also 
be considered. 

• Provide an expansion of the QS aberrations that is realistic, as per Eq. (|23|) . As there 
is no requirement that the expansion functions are orthogonal, there is considerable 
freedom allowed. However, the more terms that are required, the more difficult the 
estimation task becomes. Any method for estimation of exoplanet intensities, whether 
with long exposures or short, must do this in some form. 

• Find numerical methods for handling the quadratic system in Section 1X21 The number 
of unknowns in the linear system in Section 13.11 grows quadratically with the number 
of unknown aberration coefficients, while in the quadratic formulation, it grows only 
linearly. 
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• Provide optimal estimators of AO residual <fr r (p,t) and its error statistics, given the 
WFS data. 

• Determine the bias and variance of the elements of the matrix H and the atmospheric 
speckle function A(p,t), given the statistics of the AO residual. 

• Determine the error bars on the estimated quantities within the linear and quadratic 
formulations, given the required statistics on the necessary quantities. 

• In order to reduce memory requirements and processing times, sequential estimation 
methods in the spirit of Kalman filters should be considered. 
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Fig. 1. — (color online) Optical scheme of a coronagraphic imager. The incoming AO residual is 
</V(r, t), and the upstream static aberration is denoted <p u ( r ) (downstream static aberrations are not 
included in this analysis). The pupil plane fields upstream and downstream of the coronagraphic 
focal plane mask and downstream of the Lyot stop are given by uj(r,t), and u s (r,t), respectively. 
The aperture function A(r) accounts for both the telescope aperture and the Lyot stop. From 
Sauvage et al. (2010). 
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Fig. 2. — (color online) The black solid line shows a portion of the time-series of the temporal 
variation of the planetary intensity from the Experiment 1 in Sec. 01 The dashed blue line shows 
the corresponding stellar intensity at the planet's location in the image plane. Both the planetary 
and stellar intensity are normalized to have a mean of unity in this figure. The stellar intensity is 
enhanced by sinusoidal static aberration at the spatial frequency corresponding to planet's location 
(cf. Fig. E]). 
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Fig. 3. — (color online) Histograms of the planetary intensity (top) and the stellar intensity 
(bottom) from the full time-series excerpted in Fig. [2j These distributions exhibit the negative 
(top) and positive (bottom) skewness discussed in Gladysz et al. (2010). 
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Fig. 4. — (color online) The 2.2 /i The black curve shows Strehl ratio for the Experiment 1, and 
the red curve below shows the fraction of light removed by the ideal coronagraph. 
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Fig. 5. — The static aberration function assumed for all experiments. This function is composed of 
a linear combination of high-order Zernike polynomials plus a sinusoid corresponding to the spatial 
frequency of the planet's location. 
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Fig. 6. — The Fourier transform magnitude of static aberration function shown in Figure [5j The 
two bright spots, the lower of which is spatially coincident with the planet (at distance from the star 
of 3X/D), are due to the sinusoidal aberration. The zero-padding and other details are described 
in the text. The pixel numbers correspond to the central 61 x 61 pixels of the zero-padded array. 
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Fig. 7. — The stellar intensity averaged over all 4000 exposures in Experiment 1 without noise. The 
planet is not included and it is too faint to be seen. The colorscale is in units [photons/ms/pixel]. 
The pixel numbers correspond to the central 61 x 61 pixels of the zero-padded array. 
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Fig. 8. — (color online) Static aberrations coefficients {a^}. The red pentagrams represent the true 
values used in the simulation, resulting in Fig. The black circles represent the values estimated 
by the algorithm in Experiment 1. The error bars, with values of about 1 x 10~ 4 , are smaller than 
the plotting symbols. 
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Fig. 9. — (color online) Products of aberrations coefficients {aja^}. These quantities need not be 
estimated in quadratic formulation (Section [32]). The red pentagrams represent the true values used 
in the simulation. The black circles represent the values estimated by the algorithm in Experiment 
1. The error bars are the square-roots of the diagonal elements of the covariance of the linear 
estimator. 
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